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Abstract 

We study Semidefinite Programming, SDP, relaxations for Sensor Network Localiza- 
tion, SNL, with anchors and with noisy distance information. The main point of the 
paper is to view SNL as a (nearest) Euclidean Distance Matrix, EDM , completion 
problem and to show the advantages for using this latter, well studied model. We first 
show that the current popular SDP relaxation is equivalent to known relaxations in 
the literature for EDM completions. The existence of anchors in the problem is not 
special. The set of anchors simply corresponds to a given fixed clique for the graph of 
the EDM problem. We next propose a method of projection when a large clique or a 
dense subgraph is identified in the underlying graph. This projection reduces the size, 
and improves the stability, of the relaxation. 
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In addition, viewing the problem as an EDM completion problem yields better 
low rank approximations for the low dimensional realizations. 

And, the projection/reduction procedure can be repeated for other given cliques 
of sensors or for sets of sensors, where many distances are known. Thus, further size 
reduction can be obtained. 

Optimality/duality conditions and a primal-dual interior-exterior path following 
algorithm are derived for the SDP relaxations We discuss the relative stability and 
strength of two formulations and the corresponding algorithms that are used. In par- 
ticular, we show that the quadratic formulation arising from the SDP relaxation is 
better conditioned than the linearized form, that is used in the literature and that 
arises from applying a Schur complement. 

1 Introduction 

We study ad hoc wireless sensor networks and the sensor network localization, SNL , prob- 
lem with anchors. The anchors have fixed known locations and the sensor-sensor and sensor- 
anchor distances are known (approximately) if they are within a given (radio) range. The 
problem is to approximate the positions of all the sensors, given that we have only this 
partial information on the distances. We use semidefinite programming, SDP , relaxations 
to find approximate solutions to this problem. 

In the last few years, there has been an increased interest in the SNL problem with 
anchors. In particular, SDP relaxations have been introduced that are specific to the prob- 
lem with anchors. In this paper we emphasize that the existence of anchors is not special. 
The SNL problem with anchors can be modelled as a (nearest) Euclidean Distance Matrix, 
EDM , completion problem, a well studied problem. There is no advantage to considering 
the anchors separately to other sensors. The only property that distinguishes the anchors is 
that the corresponding set of nodes yields a clique in the graph. This results in the failure of 
the Slater constraint qualification for the SDP relaxation. We then show that we can take 
advantage of this liability. We can find the smallest face of the SDP cone that contains the 
feasible set and project the problem onto this face. 

This projection technique yields an equivalent smaller dimensional problem, where the 
Slater constraint qualification holds. Thus the problem size is reduced and the problem 
stability is improved. In addition, viewing the problem as an EDM completion leads to 
improved low rank factorizations for the low dimensional realizations. And, by treating the 
anchors this way, we show that other cliques of sensors or dense parts of the graph can 
similarly result in a reduction in the size of the problem. In addition, not treating other 
cliques this way can result in instability, due to loss of the Slater constraint qualification. 

We also derive optimality and duality conditions for the SDP relaxations. This leads 
to a primal-dual interior-exterior path following algorithm. We discuss the robustness and 
stability of two approaches. One approach is based on the quadratic constraint in matrix vari- 
ables that arises from the SDP relaxation. The other approach uses the linearized version 
that is used in the literature, and that is obtained from an application of the Schur comple- 
ment. Numerical tests comparing these two equivalent formulations of the SDP relaxation 
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are included. They show that the quadratic formulation is better conditioned and requires 
fewer iterations to reach a desired relative duality gap tolerance. These tests confirm results 
in the literature, see [Ul E], on the conditioning of the central path and the comparison of 
different barriers. 

1.1 Related Work and Applications 

The geometry of EDM has been extensively studied in the literature, e.g. [IH [TU] and more 
recently in [21 Q] and the references therein. The latter two references studied algorithms 
based on SDP formulations of the EDM completion problem. 

Several recent papers have developed algorithms for the SDP relaxation designed specif- 
ically for SNL with anchors, e.g. [SI El 13 IH [2U EJ [301 HH] - Relaxations using second order 
cones are studied in e.g. [271 EE]- 

The SDP relaxations solve a closest SDP matrix problem and generally use the l\ 
norm. The £2 norm is used in [19J, where the noise in the radio signal is assumed to come 
from a multivariate normal distribution with mean and variance-covariance matrix a 2 1, i.e. 
from a spherical normal distribution so that the least squares estimates are the maximum 
likelihood estimates. (We use the £2 norm as well in this paper. Our approach follows that 
in [2] for EDM completion without anchors.) 

Various applications for SNL are discussed in the references mentioned above. These ap- 
plications include e.g. natural habitat monitoring, earthquake detection, and weather/ current 
monitoring. 

1.2 Outline 

The formulation of the SNL problem as both a feasibility question and as a least squares 
approximation is presented in Section[21 We continue in Section[3]with background, notation, 
including information on the linear transformations and adjoints used in the model. In 
particular, since this paper emphasizes using EDM , this section provides details on distance 
geometry. In particular, we provide details on the linear mappings between EDM and 
SDP matrices. 

The SDP relaxations are presented in Section [U This section contains the details for 
the four main contributions of the paper: i.e. 

(i) the connection of SNL with EDM ; (ii) the projection technique for cliques 
and dense sets of sensors; (iii) the improved approximation scheme for locating 
the sensors from the SDP relaxation; and (iv) a numerical comparison show- 
ing the better conditioning of the quadratic formulation relative to the linear 
formulation used in the literature. 

We begin in Section l4TTl with several lemmas that describe the feasible set of the SDP relaxation, 
e.g. Lemma [4.21 provides several equivalent characterizations that show the connection with 
EDM . The key to the connection is the loss of the Slater constraint qualification (strict 
feasibility); but one can project onto the minimal face in order to obtain the Slater condition 
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and guarantee numerical stability and strong duality. This Lemma also shows the equivalent 
representations of the feasible set by a quadratic and a linear semi definite constraint. Then 
Lemma 14.31 shows that the above projection idea can be used for other cliques and dense 
subgraphs. 

The optimality and duality theory for the SDP relaxations is presented in Section We 
show that strict feasibility holds for the dual if the underlying graph for the primal problem 
is connected. 

Our primal-dual interior /exterior-point (p-d i-p) algorithm is derived in Section [HI We 
include a heuristic for obtaining a strictly feasible starting point. The algorithm uses a 
crossover technique, i.e. we use the affine scaling step without backtracking once we get a 
sufficiently large decrease in the duality gap. 

We then continue with the numerical tests in Section [7J Concluding remarks are given 
in Section [HI 



2 SNL Problem Formulation 

Let the n unknown (sensor) points be p 1 , p 2 , . . . , p n e W, r the embedding dimension; and 
let the m known (anchor) points be a 1 ,a 2 ,...,d m G M r . Let X T = [p 1 ,p 2 , . . . ,p n ], and 
A T = [a 1 , a 2 , ... , a m }. We identify a 1 with p n+l i for % = 1, . . . ,m, and sometimes treat these 
as unknowns. We now define 

P T := (p\p 2 ,...,p n ,a\a 2 ,...,a m ) = (p\p 2 , . . . , p n , p n+ \ p n+2 , . . . ,p n + m ) = (X T A T ) . 

(2.1) 

Note that we can always translate all the sensors and anchors so that the anchors are centered 
at the origin, i.e. A T <— A T — ■^A T ee T yields A T e = 0. We can then translate them all back 
at the end. In addition, we assume that there are a sufficient number of anchors so that 
the problem cannot be realized in a smaller embedding dimension. Therefore, to avoid some 
special trivial cases, we assume the following. 

Assumption 2.1 The number of sensors and anchors, and the embedding dimension satisfy 

n >> m > r, A T e = 0, and A is full column rank. 

Now define (Af e ,Af u ,Ni), respectively, to be the index sets of specified (distance values, 
upper bounds, lower bounds), respectively, of the distances dij between pairs of nodes from 
{p l }i (sensors); and let (M. e , M. u , -Mi), denote the same for distances between a node from 
{p l }i (sensor) and a node from {a fc }^ (anchor). Define (the partial Euclidean Distance 
Matrix) E with elements 

( d\ if ij e Af e U M e 

Eij = { \\P l ~P*\\ 2 = \W~ n — a j ~ n \\ 2 if i,j>n 

y otherwise. 

The underlying graph is 

G = 0>,E), (2.2) 
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with node set V = {1, . . . , m + n) and edge set £ = M e U M e U {ij : i,j > n}. Note that the 
subgraph induced by the anchors (the nodes with j > n) is complete, i.e. the set of anchors 
forms a clique in the graph. Similarly, we define the matrix of (squared) upper distance 
bounds U b and the matrix of (squared) lower distance bounds L b for ij G M u U M u and 
Mi U Mi, respectively. 

Our first formulation for finding the sensor locations p> , j < n, is the feasibility question 
for the constraints: 
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Note that the first three and the last two sets of constraints are quadratic, nonconvex, 
constraints. We added the anchor-anchor distances to emphasize that these are not special 
and can be treated in the same way as the other distances. 

The above can also be considered as a Graph Realization Problem, i.e. we are given an 
incomplete, undirected, edge-weighted simple graph. The node set is the set of of sensors 
and anchors. The weights on the edges are the squared distance between two nodes, not all 
known and possibly inaccurate. A Realization of G in W is a mapping of nodes into points 
p % in 9fJ r with squared distances given by the weights. 

Let W p , W pa , W a be weight matrices for the sensor-sensor, sensor-anchor, anchor- anchor, 
distances respectively. For example, they simply could be 0, 1 matrices to indicate when an 
exact distance is unknown or known. Or a weight could be used to verify the confidence in 
the value of the distance. The weights in W a correspond to anchor-anchor distances and are 
large, since these distances are known. In the literature, these anchor-anchor distances are 
considered as constants in the problem. We emphasize that they are equivalent to the other 
distances, and that the SNL problem is a special case of the EDM problem. If there is 
noise in the data, the exact model (I2.3P can be infeasible. Therefore, we can minimize the 
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weighted least squares error. 

min /!(/>):= \ ^ {W,)^ - p>\\ 2 - E l3 f 
(jj')eJV e 

+ 1 E (^paWlb'-^ll 2 -^ 

(j,fc)eX e 

- p>\\ 2 - 

,SNL LS ) sub j ec t to - p?'|| 2 < t/* V(i,j)eK ( ^ ^ 

y(i,k) e M u 

V(i,j)eM 

V(i,fc) eM, 
Vz, j > n) . 

This is a /iani problem to solve due to the nonconvex objective and constraints. We again 
included the anchor-anchor distances within brackets both in the objective and constraints. 
This is to emphasize that we could treat them with large weights in the objective or as 
holding exactly without error in the constraints. 





3 Distance Geometry 

The geometry for EDM has been studied in e.g. [22} [131 EEl EE], and more recently, in e.g. 
[2], p. Further theoretical properties can be found in e.g. 121 Q31 [161 GEl EDI E2[ EU] . Since 
we emphasize that the EDM theory can be used to solve the SNL , we now include an 
overview of the tools needed for EDM . In particular, we show the relationships between 
EDM and SDP . 



3.1 Linear Transformations and Adjoints Related to EDM 

(We use the notation from [TjJ]. We include it here for completeness.) We work in spaces of 
real matrices, A4 sxt , equipped with the trace inner-product (A, B) = trace A T B and induced 
Frobenius norm || = trace A T A. For a given B e S n , the space ofnxn real symmetric 
matrices, the linear transformation diag (B) e R n denotes the diagonal of B\ for v € M n , 
the adjoint linear transformation is the diagonal matrix diag*(t>) = Diag (i>) G S n . We now 
define several linear operators on S n . (A collection of linear transformations, adjoints and 
properties are given in the appendices.) 

V e (B) := diag( J B)e T + ediag( J B) T , K{B) := V e (B)-2B, (3.5) 

where e is the vector of ones. The adjoint linear operators are 

V* e (D) = 2Diag (De), K?{D) = 2(Diag (De) - D). (3.6) 
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By abuse of notation we allow T> e to act on R n : 

V e (v) = ve T + ev T , v G R n . 

The linear operator JC maps the cone of positive semidefmite matrices (denoted SDP ) onto 
the cone of Euclidean distance matrices (denoted EDM ), i.e. JC{SDP ) = EDM . This 
allows us to change problem EDMC into a SDP problem. 

We define the linear transformation sblkj(S') = Si G <S*, on S G S n , that pulls out the 
2-th diagonal block of the matrix S of dimension t. (The values of t and n can change and 
will be clear from the context.) The adjoint sblk*(T) = sBlkj(T), where T G <S*, constructs 
a symmetric matrix of suitable dimensions with all elements zero expect for the z-th diagonal 
block given by T. 

Similarly, we define the linear transformation sblkj 3 (G) = Gij, on G G S n , that pulls out 
the ij block of the matrix G of dimension k x I and multiplies it by y/2. (The values of k, 
I, and n can change and will be clear from the context.) The adjoint sblk *■(,/) = sBlky(J), 
where J G Ai kxl = M. kl , constructs a symmetric matrix that has all elements zero expect 
for the block ij that is given by J multiplied by and for the block ji that is given by 

J T multiplied by The multiplication by y/2 (or -j=) guarantees that the mapping is an 
isometry. We consider J G A4 kxl to be a k x I matrix and equivalently J G M fc/ is a vector 
of length kl with the positions known. 

3.2 Properties of Transformations 

Lemma 3.1 flT$) Define the linear operator on S n by 

offDiag (S) = S - Diag (diag (S)). 
Let J .= I — ^ee T . Then, the following holds. 

• The nullspace N (JC) equals the range lZ(V e ). 

• The range TZ(fC) equals the hollow subspace of S n , denoted Sh '■= {D G S n : diag (D) = 
0}. 

• The range TZ(JC*) equals the centered subspace of S n , denoted S c := {B G S n : Be = 0}. 

• The Moore-Penrose generalized inverse fO(D) = J (offDiag (D)) J. ■ 

Corollary 3.1 (JWj) 

1. Let So denote the subspace of diagonal matrices in S n . Then 

S c = M(D* e ) = n(tC*) = K(K)) 1 M{K) = K(V e ) 

s H = n{K)=M{v e ) l s D = M{K*) = n{vt). 
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2. Let [V 77^ e) be an nx n orthogonal matrix. Then 

YhO Y = VYV T + V e (v) h 0, for some Y e S"-\v E W 1 . 



Let B = PP T . Then 



D i:j = Wrf-pif = (diag( J B)e T + ediag {B) T - 2B) i5 = {K{B)) 



(3.7) 



i.e. the EDM D = {D^) and the points Pi in P are related by D = JC(B), see (13.51) . 



4 SDP Relaxations of SNL based on EDM Model 



We first study the SDP relaxation used in the recent series of papers on SNL , e.g. [6], HI 
[2HE1IT7]. (See (14.91) and Section |4\ 1 . 31 below. ) This relaxation starts by treating the anchors 
distinct from the sensors. We use a different derivation and model the problem based on 
classical EDM theory, and show its equivalence with the current SDP relaxation. By 
viewing the SNL problem as an EDM problem, we obtain several interesting results, e.g. 
clique reduction, and a geometric interpretation on how to estimate sensor positions from 
the SDP relaxation optimum. 

4.1 Connections from Current SDP Relaxation to EDM 

Let Y = XX T . Then the current SDP relaxation for the feasibility problem for SNL uses 



Lemma 3.2 (H3j) Suppose that ^ B G <S n . Then D = K(B) is EDM. 




(4.8) 



This is in combination with the constraints 
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(4.9) 



trace 



Z, 



Eij, Vzj G A4 e , i < j = n + k. 
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4.1.1 Reformulation using Matrices 

/ XX T XA T \ 

We use the matrix lifting or linearization Y := PP T = ( j^-^t j^j^t ) an d Z := [I; P] [I; P] T 
I P T \ 

p y j ' ^ e dimensions are: 

X e M n ' r ; A e M m ' r ; P e M m+n ' r ; Y e S m+n - Z e s rn+n+r . 

Adding the hard quadratic constraint Y = PP T allows us to replace the quartic objective 
function in SNLls with a quadratic function. We can now reformulate SNL using matrix 
notation to get the equivalent EDM problem 

min f 2 (Y):=±\\Wo(lC(Y)-E)f F 

subject to g u (Y) :=H u o (JC(Y) - U b ) < 

(SNL M ) g l {Y):=H l o{K,{Y)-L b ) > (4.10) 

Y-PP T = 

(sblk 2 (/C(y)) = IC(AA T )) , 

where W G S n+m is the weight matrix having a positive ij-element if G f\f e UM e U{{ij) : 
i,j > n}, otherwise. H u , Hi are 0, 1-matrices where the ij-th element equals 1 if an upper 
(resp. lower) bound exists; and it is otherwise. By abuse of notation, we consider the 
functions g u ,9i as implicitly acting on only the nonzero components in the upper triangular 
parts of the matrices that result from the Hadamard products with H u ,Hi, respectively. 
We include in brackets the constraints on the clique corresponding to the anchor-anchor 
distances. 

Remark 4.1 The function f2(Y) = f 2 (PP T ), and it is clear that f 2 (PP T ) = fi(P) in 
(12. 4p . Note that the functions f2,g u ,9i act only on Y and the locations of the anchors 
and sensors is completely hiding in the hard, nonconvex quadratic constraint Y = PP T = 
( XX T XA T \ 

( j^x T AA T ) ' ^ 6 "P ro ^ em SNLm is a linear least squares problem with nonlinear 

constraints. The objective function is generally underdetermined. This can result in ill- 
conditioning problems, e.g. [77] /. Therefore, reducing the number of variables helps with 
stability. 

4.1.2 SDP Relaxation of the Hard Quadratic Constraint 

We now consider the hard quadratic constraint in (I4.10p 



where P is defined in ( 12. ip . We study the standard current semidefinite relaxation in (14.91) 
with (14. 8p . or equivalently with Y >z PP T ■ We show that this is equivalent to the simpler 
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Y >z 0. We include details on problems and weaknesses with the relaxation. We first present 
several lemmas. We start with the following well known result. We include a proof for 
completeness. 

(Y Y t \ 

Lemma 4.1 Suppose that the partitioned symmetric matrix ( AA T ) — ^' Then ^21 = 
XA T ; with X = YlA(A T A)- 1 . 

Proof. Let A = WE r V T be the compact singular value decomposition, ~< H r E S r . And, 
suppose that (U U ) is an orthogonal matrix. Therefore, the range spaces 7Z(U) = 71(A) 
and the nullspace M(A T ) = 1Z(U). Consider the nonsingular congruence 



~< 




( z 


Y T 
1 21 


\Y 21 


AA T 


Yi(U 


U) 


(1 


I) 







I 

(U u 



This implies that Y^U = 0. This in turn means that MiY^) D Af(A T ), or equivalently, 
71(A) D 7l(Y 21 ). Note that the orthogonal projection onto 71(A) is A(A T A)- 1 A T . Therefore, 
Yl = Y 2 T 1 A(A T A)~ l A T = (Y 2 \A(A T A)- 1 ) A T } i.e. we can choose X = Y^_A(A T A)~ 1 . ■ 



In the recent literature, e.g. [TJ El [T7j, it is common practice to relax the hard constraint 
(14. lip to a tractable semidefinite constraint, Y >z PP T , or equivalently, Yu >z XX T with 
Y21 = AX T . The following lemma presents several characterizations for the resulting feasible 
set. 

Lemma 4.2 Let A = UT, r V T be the compact singular value decomposition of A, and let 
P, Y be partitioned as in (j2.1j) . (j4.1ip . 

\Pj \Y 21 Y 22 J- 

Define the semidefinite relaxation of the hard quadratic constraint (14.111) as: 

G(P, Y) := PP T -Y^O, Y 22 = AA T , P 2 = A. (4.12) 

By abuse of notation, we allow G to act on spaces of different dimensions. Then we get the 
following equivalent representations of the corresponding feasible setTc- 

T G = {(P, Y) : G(P, Y) z< 0, Y 22 = AA T , P 2 = A} (4.12a) 



T G = { (P,Y) : G(X,Y) * 0,y u = Y,Y 21 = AX , Y 22 = AA , P 



(4.12b) 



11 



J, X — Z 21 , P 



J 22 — 1 ) ^ — ^21 



.4 



(4.12c) 



= J } (P. V) : r b 0. 1',, = AA T , X — Y£A(A T A)-\P= [ X a )\ ( 1.12(1) 




r<,= {(P.v):z= [ Z 7 U z 7 21 |Mi.r 

^21 ^22 
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l Z (o 
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S?, X = Y£A(A T A)-' = Zl^V T , P=( 



(4.12e) 



Moreover, the function G is convex in the Lowner (semidefinite) partial order; and the 
feasible set JF G is a closed convex set. 

Proof. Recall that the cone of positive semidefinite matrices is self-polar. Let Q y and 
4>q{P) = trace Q PP T . Convexity of G follows from positive semidefiniteness of the Hessian 
V 2 4>q(P) = I <E> Q, where ® denotes the Kronecker product. 
In addition, 

oy G (P,Y) = pp T -Y=( x A f T Zy L 2 1 i XAT ~ Yi ) 

holds if and only if 

t G(X, Y n ) = XX T - Y u , and AX T - Y 21 = 0. 

This shows the equivalence with ( I4.12bl) . A Schur complement argument, with Yu = Y, 
shows the equivalence with ( X J y 0, i.e. with the set in f l4.12cl) . The equivalence 



with (14. 12dl) follows from Lemma [4.11 

To show the equivalence with the final expression (14.12el) . we note that Y y 0, Y22 = AA T , 
implies that there is no strictly feasible Y y 0. Therefore, we project the feasible set onto 
the minimal cone or face (see [8]). This yields the minimal face that contains the feasible 
set of Y, i.e. 

The result follows since the constraint P22 = AA T holds if and only if Z w is blocked as 

^ j-,2 ) — 0- (More simply, one can show the equivalence of (14. 12e[) with f )4.12cl) 

by using the compact singular value decomposition of A. However, the longer proof given 
above emphasizes that the reduction comes from using a projection to obtain the Slater 
constraint qualification.) ■ 
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The above Lemma 14.21 shows that we can treat the set of anchors as a set of sensors 
for which all the distances are known, i.e. the set of corresponding nodes is a clique. The 
fact that we have a clique and the diagonal m x m block AA T in Y is rank deficient, r < m, 
means that the Slater constraint qualification, Y y 0, cannot hold. Therefore, we can project 
onto the minimal cone containing the feasible set and thus reduce the size of the problem, 
see Lemma 14.21 (14.12e[) . i.e. the variable Y G S n+m is reduced in size to Z 6 S n+r . The 
reduction can be done by using any point in the relative interior of the minimal cone, e.g. 
any feasible point of maximum rank. The equivalent representations in (14.12cj) and in (I4.12el) 
illustrate this. 



4.1.3 Current SDP Relaxation using Projection onto Minimal Cone 

The above reduction to Y in Lemma 14.21 (14.12bl) , allows us to use the smaller dimensional 
semidefinite constrained variable 

' A ' * y e S n+m , Y u = Y, Y 21 = AX T . (4.19) 



X Y 

This is what is introduced in e.g. [6]. 

Remark 4.2 Note that the mapping Z s = Z S (X,Y) : M n,r x,S"-> S n+r is not onto. This 
means that the Jacobian of the optimality conditions cannot be full rank, i.e. this formulation 
introduces instability into the model. A minor modification corrects this, i. e. the I constraint 
is added explicitly. 

r/ 1 J 21 ]h0, Z n = I,Yu = Z 2 2,Y 2 \ = AZ^x- 

Zj21 Z-22 / 

To develop the model for computations, we introduce the following notation. 
x := vec ^sblk 2 i ^ =V / 2vec(X), y:=svec(Y), 

where we add a/2 to the definition of x since X appears together with X T in Z s and implicitly 
in Y, with Y 2 \ = AX T . We define the following matrices and linear transformations: 

Z x s (x) : = sBlk 2i (Mat (x) ) , Z y s (y) : = sBlk 2 (sMat (y) ) , 
Z s (x,y) := Z*{x) + Zy(y), Z s := sBlk ^7) + Z s (x, y), 

y x (x) := sBlk 21 (AMat (x) T ) } _y y (y) := sBlk^sMat (y)) 
y{x, y) := y x {x) + yy{y), Y := sBlk 2 {AA T ) + y{x, y). 

E ■= Wo [E- /C(sBlk 2 (AA T ))] , 
U b := H u o [/C(sBlk 2 (Ayi r )) - U b ] , 
L b :=H l o [L b - /C(sBlk 2 (A4 T ))] . 
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By abuse of notation, we let the functions . . . , y act directly on the matrices X, Y. The 
meaning will be clear from the context. 

The unknown matrix Y in (I4.10p is equal to y(x, y) with the additional constant in the 
2,2 block, i.e. our unknowns are the vectors x,y which are used to build Y and Z s . Using 
this notation we can introduce the following vector form of the relaxation of (|4.10p . 

min f 3 (x, y):=l\\Wo (K(y(x, y))) - E\\} 

(qnt \ subject to g u (x,y) := H u o IC(y(x,y)) - U b < , 

{SNLMV) gi (x,y):=L»- Hl oK:(y( x ,y)) < ^ 

sBlk^+Z^y) t 0. 

As above, we consider the functions g u ,9i as implicitly acting only on the nonzero parts 
of the upper triangular part of the matrix that results from the Hadamard products with 
H u ,Hi, respectively. 

4.1.4 SDP Formulation Using EDM 

The equivalent representations of the feasible set given in Lemma |4~2| in particular by (14. 12ej) . 
show that SNL is an EDM problem D = JC(Y), with the additional upper and lower 
bound constraints as well as the block constraint sBlk 2 (.D) = JC(AA T ), or equivalently, 
sBlk 2 CP) = A4 T . 

Remark 4.3 Suppose that we can increase the size of the clique containing the anchor nodes 
by adding sensor nodes where the distances are exactly known. Then these sensor nodes can 
be treated as anchor nodes, though their position is unknown. 

We can now obtain an equivalent relaxation for SNL by using the EDM (14.101) and 
replacing the hard quadratic constraint with the simpler semidefmite constraint Y >z 0. We 
then observe that the Slater constraint qualification (strict feasibility) fails. Therefore, we 
can project onto the minimal cone, i.e. onto the minimal face of the SDP cone that contains 
the feasible set. see 012]. Let 

Ua= ( I a)' Us= ( I u) ■' Where A = U ^V T . (4.21) 

Recall that sblk 2 (/C(f/ s Z[/J)) = JC(AA T ) is equivalent to Z 2 2 = S^. We get two SDP relaxation 
that are equivalent to SNLmv '■ 

min U{Z) := \\\W o (JC(U s ZUj) - E)\\ 2 p 
subject to H u o (IC(U s ZUj) - U h ) < 



[SNL ED ms) H t o (K,{U S ZUJ) - L b ) > (4.22) 

Z22 = 2;? 

z y 0. 
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and 

min f 3 (Z) := ||| W o {K,{U A ZUl) - E)f F 
subject to H u o (K.(U A ZUl) - U b ) < 

(SNL EDM A ) Hi o (IC(U A ZUl) - L b ) > (4.23) 

Z22 = 

z y 0. 



Remark 4.4 iVote that we do not substitute the constraint on Z 2 2 into Z , but leave it explicit. 
Though this does not change the feasible set, it does change the stability and the dual. This 
can be compared to the SDP relaxation for the Max-Cut problem with constraint that the 
diagonal of X is all ones, diagX = e and X y 0. However, one does not substitute for the 
diagonal and rewrite the semidefinite constraint. 



4.2 Clique Reductions using Minimal Cone Projection 

Now suppose that we have another clique of p > r sensors where the exact distances are 
known and are used as constraints. Then there exists a matrix Y = PP T that has a diagonal 
rank deficient p x p block. Since all feasible points are found from elements in the set 
Y + A/"(/C), we conclude that for p large enough, the diagonal block remains rank deficient 
for all feasible Y, i.e. the Slater constraint qualification fails again, if the corresponding 
distances are added as constraints. 

We now see that we can again take advantage of the loss of the Slater constraint qualifi- 
cation. 



Lemma 4.3 Suppose that the hypotheses and definitions from Lemma \4^ hold; and suppose 



that there exists a set of sensors, without loss of generality S c := {p t+1 , . . . ,p n }, so that the 
distances \\p l — p*\\ are known for all t + 1 < i,j < n; a k , i.e. the graph of the partial 
EDM has two cliques, one clique corresponding to the set of known anchors, and the other 
to the set of sensors S c . Let P, Y be partitioned as 

P = I Pa I = I . I , Y = Kn %>. YZ I = PP 1 




A 




132 
^33 



where Pi = Ai,i = 2,3, and A 3 = A corresponds to the known anchors while P2 = A 

P 

EDM , E = K{Y), be correspondingly blocked 



corresponds to the clique of sensors and X — [ ^} \ corresponds to all the sensors. Let the 



E 1 ■ 

E = I ■ En 



E, 
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so that E 3 = K(AA T ) are the anchor- anchor squared distances, and E 2 = )C(P 2 P 2 ) are 
squared distances between the sensors in the set S c . Let 

B = K)(E 2 ). 

Then the following hold. 
1. Be = and 

Y 22 = B + y 2 e T + ey 2 h 0, for some y 2 G 71(B) + ae, a > 0, with rank (Y 22 ) < r. 

(4.24) 



2. The feasible set Tq in Lemma \4 . 2\ can be formulated as 



Tc; := { ( P. Y) Z = \ Z 2J Zoo 4 I [ U 2 |Z|0 ( r 2 





fZ n 


7 T 
z 21 


7 T 
^31 


H 


Z21 


Z22 


7 T 
^32 






^32 


■^33 




Z33 




X = 




or equivalently as 



(4.25) 



F G ={(P,Y):Z=\ Z 21 Z 22 Z 3 T 2 |^0,F=|0 *7 2 0|Z|0 C/ 2 

yT 
^31 

C/ 2 Z 32 

(4.26) 

where B := B + 2ee T = (U 2 U 2 ) ( ^ 2 ^ ) ( t/2 U 2 ) T is the orthogonal diagonal- 







7 T 

Zj 2\ 


2§\ 


H 


Z 2 \ 


Z22 


-^32 




\z 31 


Z« 


-^33 / 






■^33 = 


L r , X 






ization of B, with D 2 G S+ ,r 2 < r + 1. 



Proof. We proceed just as we did in Lemma I4T2| i.e. we reduce the problem by projecting 
onto a smaller face in order to obtain the Slater constraint qualification. 

The equation for Y 22 for some y 2 , given in (14.241) . follows from the nullspace character- 
ization in Lemma I3.ll Moreover, Y 22 = P2P 2 implies that rank (Y22) < r, the embedding 
dimension. And, Y 22 >z 0, Be = implies the inclusion y 2 G 71(B) + ae, a > 0. Moreover, 
we can shift P 2 = P 2 — ^—^(P 2 e)e T . Then for B = P 2 P 2 , we get Be = 0, i.e. this satisfies 
B = lO(E 2 ) and rank (B) < r. Therefore, for any Y = B + ye T + ey T >z 0, we must have 
y = ae, a > 0. Therefore, B has the maximum rank, at most r + 1, among all feasible 
matrices of the form __< Y G B + M(K) . B determines the smallest face containing all such 
feasible Y. 

Define the linear transformation H : M n ~' — > 5 n_t by H(y) = Y 22 + ye T + ey T . Let 
C := Y 22 + 7Z(7) e ) and jF e denote the smallest face of iS™~* that contains C Since I? is a 
feasible point of maximum rank, we get 

B = B + V e (y 2 ) G (C n relint .F e ) . 
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Thus the face 



T e = {U 2 ZUZ : Z E S r + 2 } = {Y E 5+"* : trace Y{U 2 U%) = 0}. 
Now, we expand 

Yn Yji Y_l\ (h 

Y 2 \ Y 22 Y^ J = 10 U 2 
Y m Yo K, / V 




U 2 Z 2 i U 2 Z 22 U2 U 2 Z 32 V S 
U S Z 31 V S Z 32 U^ U s Z 2 r UT 



Therefore, — I %tttt )• Therefore, the expression for Z 33 and X in (I4.25P follows 

\U 2 Z 32 U S J 

from equation f)4.12el) in Lemma 14.21 The result in (I4.25P can be obtained similarly or by 
using the compact singular value decomposition of A. ■ 



Remark 4.5 The above Lemma \4.3\ can be extended to sets of sensors that are not cliques, 
but have many known edges. The key idea is to be able to use (Wi o K)^(Wi o Ej) and to 
characterize the nullspace ofWi o JC. This is studied in a forthcoming paper. 

We can apply Lemma I4T31 to further reduce the SDP relaxation. Suppose there are a group 
of sensors for which pairwise distances are all known. This should be a common occurrence, 
since distances between sensors within radio range are all known. Without loss of generality, 
we assume the of sensors to be {p t+1 , . . . ,p n }. Let E 2 , B = K){E 2 ), and U 2 , be found using 
Lemma [4.31 and denote 



Via 



(In 









(In 





;: 




v 2 





, V 2s := | 









\ o 





A) 









v 



(4.27) 



In SNL 

edm si we can replace V s with U 2s and reach a reduced SDP formulation. Simi- 
larly, for SNL 

edm a- Furthermore, we may generalize to the k clique cases for any positive 
integer k. We similarly define each Ui, 2 < i < k, and define 



VkA 



(In 




\o 





v 2 














v k 








A) 



Vu 



(In 

v 2 












... \ 

... 

V k 

V ) 



(4.28) 



17 



Then we can formulate a reduced SDP for k cliques: 



min U{Z) := §|| W o (IC(U kA ZUl A ) - E)\\ F 
subject to H u o (K,(U kA ZUl A ) - U b ) < 

(SNL k 

—cliques— A ) H t o {1C{U kA ZUl A ) - L b ) > (4.29) 

where Z kk is the last r by r diagonal block of Z. Similarly, we get 

min U{Z) :=\\\Wo(lC{U ks ZUl)-E)\\ 2 F 
subject to H u o (JC(U ks ZU£ s ) -U b )<0 

(SNL k - cUques - s ) E { o (lC(U ks ZUl) - L b ) > (4.30) 

Z kk = I r 

z^o 

For a clique with r e sensors, a f/j is constructed with r e rows and at most r + 1 columns. 
This implies the dimension of Z has been reduced by r e — r — 1. So if r = 2, cliques 
larger than a triangle help reduce the dimension of Z. As mentioned above, the existence of 
cliques is highly likely, since edges in the graph exist when sensors are within radio range. 
Moreover, the above technique extends to dense sets, rather than cliques. The key is finding 
B = (W o JCY(W o En), for an appropriate submatrix En, as well as deriving the nullspace 
ofWolC. 

4.3 Estimating Sensor Positions based on EDM Model 

After we solve SNL EDM A (or equivalently SNL EDMs ) to get an optimal solution Z s , we 

can express 



Y = U s Z s Uj = ( X) 1 Y J X ) , ^12 = AA T , Y 21 = AX T , for some X. 



21 

Y2I Y 2 2 

To complete the SNL problem, we have to find an approximation to the matrix P G 

Ji4 n+m,r , i.e. the matrix that has the sensor locations in the first n rows, also denoted 

X, and the anchor locations in the last m rows, denoted A. 

1 P P 

m+m ™ n +r, E> _ / . n . 21 



Since Y G Sl +m , there exists P = 5 ~ 21 G M n+m such that PP T = Y. By 

\ P. 12 r 22 / 

Assumption 12.11 the anchors are centered, i.e. A T e = 0. We can translate the locations in 
P, so that the last m locations are centered, i.e. without loss of generality we have 

(A2 A 2 ) T e = 0, PP T = Y. (4.31) 

Also 

{P G M n+m : Y = PP T ) = {PG M n+m : P = PQ, for some orthogonal Q G M n+m } . 
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In other words, from the optimal Y, all the possible locations can be obtained by a rota- 
tion/reflection of P. However, these locations in the rows of P are in M n+m , rather than in 
the desired embedding space M r , where the anchors lie. 

Remark 4.6 Since SNL is underdetermined, in general, the optimum Y is not unique. 
Therefore, finding a lower rank optimum Y should result in better approximations for the 
sensor locations. 

Following are two methods for finding an estimate to the sensor locations, X. The first 
is the one currently used in the literature. The second is a strengthened new method based 
on the EDM interpretation. 

In the recent papers on SNL, e.g. [7J [6j [17], X is taken directly from the optimal 

see e.g. (I4.19p . Equivalent ly, since A is full column rank r, and 



X Y 

the equations in AX T = Y 2 \ are consistent, we can solve for X uniquely from the 
AX T = Y21. We now describe the underlying geometry of using this X. 

Recall that A = UE r V T and (P u P 22 ) ( A 2 A2 f = AA T = (A 0)(A 0) T . 
Therefore, these three matrices all have the same spectral decomposition and all can be 
diagonalized using U. This implies that the three matrices ( P12 P22 ) > A, ( A ) can 
all use the same set of left singular vectors in a compact singular value decomposition, 
SVD. Therefore, ( A2 P22 ) Q = (A ), for some orthogonal Q, i.e. 

3Q, Q T Q = I, with P = PQ = ( P £ (4.32) 

This yields 

Pll Pl2\ y _ ppT _ ( Yn P11A 



Since AP^ = Y21 = AX T , we see that Pn = X. Thus the first n rows of P project 
exactly onto the rows of X, after the rotation/reflection with Q to make the bottom m 
rows equal to A. If we denote the orthogonal projection onto the first r coordinates by 
P r , then the resulting operation on the locations in the rows of P can be summarized 
by 

P r p T = (p r Q T \ P T G W ® {0} n+m ~ r , Y w Y p := PP r P T . 

Note that the product P r Q T is not necessarily idempotent or symmetric, i.e. not 
necessarily an (orthogonal) projection. Moreover, the term that is deleted P12 can be 
arbitrary large, while the rank of Y can be as small as r + 1. The relaxation from 
Y n = XX T to Y u = ( P n P 12 ) ( P u P 12 f = AiPfi + t XX T } shows that 

using X = Pn has an error of the order of ||Pi2|| 2 . 

Method 1: Estimate the location of the sensors using X in the optimal Z s 
or, equivalently, solve for X using the equation AX T = Y21, where Y21 is 
from the optimal Y. 
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2. In Method 1, the matrix PP r P T provides a rankr approximation to Y. However, if 
II A2II i n (14.331) is large, then it appears that we have lost information. It is desirable to 
keep as much of the information from the high dimensional locations in P as we can, 
i.e. the information that is contained in Y\\. If we do not consider the anchors distinct 
from the sensors, then we would like to rotate and then project all the rows of P onto 
a subspace of dimension r, i.e. we consider the problem to be an EDM completion 
problem and would like to extract a good approximation of the positions of all the 
nodes. Since the last m rows corresponding to the anchors originated from a clique, 
the corresponding graph is rigid and the corresponding projected points will be close 
to the original anchor positions. We realize this using the spectral decomposition. (See 
e.g. [2], where error estimates are included.) 

Y=(U r Ur)(^ v Wr % f . 
\ U ^-"n+m—r J 

Then, considering the problem as an EDM completion problem, we first find a best 
rank r approximation to Y, denoted Y r := U r Y> r Uj ' . Only then do we find a particular 
full rank factorization P r G M n+m ' r such that Y r = P r P?, i.e. P r = U r T}J 2 . It 
remains to find an orthogonal Q in order to find P = P r Q- Fortunately, we can use 
the information from the anchors to find the orthogonal Q. 

Method 2: Suppose that P r = U r Y> r is found as above, with P r = 

We find Q as a minimum for mingTQ = / \\P2Q — A\\ 2 F . The solution is given 
analytically by Q = VqU%, where UqUqVq = A T P 2 is the SVD for A T P 2 . 
Then the rows of P±Q are used to estimate the locations of the sensors. 

Numerical tests for the two methods, are given in Section 17.11 Method 2 proved to be 
consistently more accurate. However, method 1 locates all sets of sensors that are uniquely 
localizable in R r , see |24j . 

Remark 4.7 As above, suppose that Y is an optimum for the SDP relaxation. The problem 
of finding a best P to estimate the sensor locations is equivalent to finding 

P* e argmin , \\W o (JC(PP T - Y)) \\ F . 

Ha) 

Equivalently , we want to find 

Y* G Sl +m , rank (Y*) = r, Y* 2 = AA T , Y* = Y + N{W o /C). 

However, finding such a Y* is equivalent to finding the minimal rank matrix in the inter- 
section of the semidefinite cone and an affine space. This is still an open/hard problem. 
Recently, I23\ l21f proposed randomization methods for SDP rank reduction. These methods 
can generate a low rank positive semidefinite matrix in an approximate affine space. 
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5 Duality for SNL with Quadratic Constraint 

Instead of using the standard linearized relaxation as SNLmv in (I4.20p and in [19J, we now 
study the new relaxation without linearizing the quadratic constraint XX T — Y ^ 0. This 
avoids ill-conditioning caused by this linearization, see Remark 14.21 Our numerical results 
indicate that the new quadratic approach is more stable than the linear approach, see more 
in Section [7J A discussion on the strengths of the corresponding barriers is given in [T5J [5]. 
Recall that x = ^/2vec (X), y := svec (Y). We begin with the reduced problem 



mm 



f 3 (x,y) :=MWo(lC(y(x,y)))-E\ 



Pll 2 

2||i-f ~ v^v^ K^i n ) ) ) \\F 

/ c< t\t t \ subject to g u (x,y) := H u o IC(y(x,y)) - U < . , 

g l ^ V )^L-H l oK{y{x^)) < (5 ' 34) 



|Mat (x)Mat (x) T - sMat (y) z< 0. 



Then the Lagrangian is 

L(x,y, A u ,A h A) = ||| W o )C(y(x, y)) - E\ \ F + (A u , H u o JC(y(x, y)) - U) 

+ (A l ,L-H l o)C(y(x,y))) (5.35) 

+ (A, ±Mat (z)Mat (x) T - sMat (y)\ , 



where < A u , < A/ G S m+n , and ^ A G S n . In addition, we denote 

X u : = svec (A u ), Xi := svec (A;), 
h u := svec (H u ), hi : = svec (Hi), A := svec (A). 

And, for numerical implementation, we define the linear transformations 

h™ = svec U (H U ) G R nz \ h? z = svec ,(#,) G R nZl , (5.36) 

where h™ z is obtained from h u by removing the zeros; thus, nz u is the number of nonzeros 
in the upper-triangular part of H u . Thus the indices are fixed from the given matrix H u . 
Similarly, for hf z with indices fixed from Hi. We then get the vectors 

X™ = svec u (A u ) G R nZu , X? z = Bvec,(A,) G R nzi . 

The adjoints are sMat u , sMat and, for any matrix M we get 

H u oM = sMat u svec U (H U o M). 

This holds similarly for Hi o M. Therefore, we could rewrite the Lagrangian as 

L(x,y,A u ,A h A) = L(x,y, A™ 2 , A™ z , A) 

= l\\W o )C(y(x, y)) - Ef F + <svec u (A M ), svec u (#« o /C(^(x, y)) - U)) 
+ (svec z(Ai), svec , (£ - if, o /C(^(x, y)))> 

+ ^A, |Mat (a;) Mat (x) T - sMat (y)^ . 

(5.37) 
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To simplify the dual of SNLmn , i-e. the max-min of the Lagrangian, we now find the 
stationarity conditions of the inner minimization problem, i.e. we take the derivatives of L 
with respect to x and y. We get 

= V x L(x,y,K u ,K h K) 

= [W o (icy*)]* (W o K(y(x, y)) -E) + [H u o (Ky*)Y (A u ) (5.38) 

- [H t o (JCy x )}* (A,) + vec (AMat (x)). 

Note that 



T(x) = ^A, | Mat (a;) Mat (x) , ( . OQ) 
= \ {x, vec (AMat (x))) . 

Therefore, dT }^ = vec (AMat (x)), since (x, vec (AMat (#))) is a quadratic form in x. Simi- 
larly, 

= V y L(x, y, A u , A/, A) 

= o (/C^)]* o /C(3^(x, y)) - £) + [H u o (/Cyv)]* (A u ) (5.40) 

- [H t o (/C^)]* (A/) - svec (A), 

since (A, sMat (y)) is linear in y. We can solve for A and then use this to eliminate A in the 
other optimality conditions, i.e. we eliminate t(n) variables and equations using 

svec (A) = [Wo(Kyv)]*(WoK(y(x,y))-E) + [H u o(1Cyv)]*{K) ( 

-[^o(/C^)]*(A z ). l ^ 4iJ 

We now substitute for A in the first stationarity condition (j5.38p . i.e. 

= [f o (^F)]* o + [5 U o (KF)r (i) - Pi ° (^)F (Ai) 

+vec (sMat { [W o (Ky v ))* (W o JC{y{x, y)) - E) 
+ [H u o {Kyv)\* (A w ) - o (Kyy))* (A,)} Mat (x)) . 

(5-42) 

The Wolfe dual is obtained from applying the stationarity conditions to the inner min- 
imization of the Lagrangian dual (max-min of the Lagrangian), i.e. we get the (dual 
SNL MN ) problem 



max L(x,y , \ u ,\ h \) 
(JE38]),(J52QD 
sMat (A u ) > 
sMat (A) >z 0. 



'IKTT n\ sub j ectto ^M,^M 

y SNL MV - D ) gMat > Qj ^ > Q (5.43) 



We denote the slack variables 

S, 

z 



U - H u o (K. (y(x, y))) , s u = svec S u 

H l0 (JC(y{x,y)))-L, Si = svec Si (5.44) 

y - xx T y o. 



We can now present the primal-dual characterization of optimality. 
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Theorem 5.1 The primal-dual variables x, y, A, X u , A/ are optimal for SNLmn if an d only 
if: 

1. Primal Feasibility: 

s u > 0, si > 0, in (EBD, 

1 T 

-Mat (x)Mat (x) - sMat (y) ^ 0. (5.45) 

2. Dual Feasibility: Stationarity equations (15.381) . ( 15. 40j) hold and 

A = sMat (A) y 0; A u > 0; A, > 0. (5.46) 

3. Complementary Slackness: 

An O S u = 

A ; o Si = (5.47) 

AZ = 0. 



We can use the structure of the optimality conditions to eliminate some of the linear dual 
equations and obtain a characterization of optimality based mainly on a bilinear equation 
and nonnegativity/ semidefmiteness. 

Corollary 5.1 The dual linear equality constraints (15.401) in Theorem \5.1\ can be eliminated 
after using it to substitute for A in (15.381) . i.e. we get equation ( 15.42ft . The complemen- 
tarity conditions in (15.471) now yield a bilinear system of equations F(x, y, X u , A/) = 0, with 
nonnegativity and semidefinite conditions that characterize optimality of SNLmn ■ M 



6 A Robust Primal-Dual Interior-Point Method 

We now present a primal-dual interior-point method for SNLmn , see in [19] for the lin- 
earized case, SNLmv ■ First, we define the equation (15.421) to be: 

L s (x,y, An, Aj) = 0. 

Then, to solve SNL MN we use the Gauss-Newton method on the perturbed complementary 
slackness conditions (written with the block vector notation): 



FJx,y,X u ,Xi 



( An o s u - fi u e ^ 
A; o si - fiie 
AZ - fi c I 

\ L s ) 



(6.48) 
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where s u = s u (x,y), s t = s t (x,y), A = A(x,y, X u , A/), Z = Z(x,y) and L s = L s (x, y, A u , A/). 
This is an overdetermined system with 

(m u + n u ) + {mi + n,[) + n 2 + nr equations; nr + t{n) + {m u + n u ) + {mi + ni) variables. 
6.1 Linearization 

We denote the Gauss-Newton search direction for (16.481) by 

f Ax \ 



The linearized system for the search direction As is: 

F'^As) s F'^x, y, X u , A,) (As) = -F^x, y, X u , A,). 

To further simplify notation, we use the following composition of linear transformations. 
Let H be symmetric. Then 

K* H (x) := ffo(K(f(i))), 
/C£(v) := if o(/C(;V(y))), 
M*,2/) := o(/C(y(x,y))). 

so, we have the following: 

A(x, y, X u , Xi) = sMat [{K? w )*{K w {x, y) - E) + (/C^)*(sMat (X u )) - (/C^)*(sMat (A,))], 
L s (x,y, A„,A Z ) = (/C^)*(i%(z,y) - £) + (K&J * (sMat (X u )) - (K%)*(sMat (A*)) 

+vec (AMat (x)). 

Define the linearization of above functions as: 

AA(Ax, Ay, AX U , AXi) = sMat [{K v w )*{K w {Ax, Ay)) + (/C^)*(sMat (AA„)) 

-(X&,)*(sMat (AA,))], 

AL s (Ax, Ay, AA U , AA;) = (/C^)*(K w (Ax, Ay)) + (/CfJ*(sMat (AA U )) - (/C^)*(sMat (AA,)) 

+vec (A AMat (x)) + vec (AMat (Ax)). 

The linearization of the complementary slackness conditions results in four blocks of equa- 
tions 
1. 

-A u o svec fC Hu (Ax, Ay) + s u o AX U = a u e - X u o s u 

2. 

Xi o svec/Ctf ; (Ax, Ay) + s/ o AA; = \h X e - X\ o s t 
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3. 



4. 



and hence 



A(sMat (Ay) - ^Mat (x)Mat (Ax) T - ^Mat (Ax) Mat (x) T ) 

Z Z 

+AA(As)(sMat (y) - -Mat (x)Mat {x) T ) 

Z 

H C I - AZ 



AL s (As) = -L s (x,y, A u ,A t 



( -\ u o svec JC Hu (Ax, Ay) + s u o AA U \ 
A; o svec JCjjt (Ax, Ay) + si o AXi 

A(sMat (Ay) - ±Mat (x)Mat (Ax) T - ±Mat (Ax) Mat (x) T ) 
+AA(As)Z 

\ AL s (As) ) 



where F' : M nxr x x sft*( m +") x sft*( m +") -> sft*(™+") x K*( m +") x _M nx ™ x K nr , i.e. the 
linear system is overdetermined. 

We need to calculate the adjoint (F^)*. We first find (1C X H )*, ()C y H )*, and ()C H )*. By the 
expression of y(Ax, Ay), we get 



(y x )*(S) \ _ ^Mat*(sblk 21 (S) T A) 



y*(S) y(yy)*(s) J V sMat*(sblki(5)) 
By the expression of ICh(Ax, Ay), we get 



vec(sblk 21 (S) T A) 
svec (sblk i(S)) 



(6.49) 



^h(^) - 1 (^)*(5) 



(F)*(/C*(i/oS)) 

(yy)*(ic*(H o s)) 



(6.50) 



Moreover, 



(AsMat(Ay), W 3 ) = trace (W^A) sMat (Ay) 



-svec (AW 3 + W 3 T A),Ay 
z 



Similarly 



-AMat(x)Mat(Ax) T , W3 ) = trace -W 3 T AMat (x) Mat (Ax 



-vcc (W / 3 T AMat(x)),Ax 
z 
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and 



-AMat(Aa;)Mat(x) T ,W / 3 



trace - W 3 T AMat (Ax) Mat (x) 



vec (AW 3 Mat (x)),Ax 



We also need to find (AA)*(S) by the expression of ICh(Ax, Ay), where S G S n , we get 

/ (K,* w )*(K, y w (svec(S))) \ I Ax \ 

(AA)*(5) = 



(/C^)*(/C^(svec(5))) 

svec[(/C»J(svec(S))] 
\-svec[(lC y Hi )(svec(S))}J 



Ay 
AX V 



(6.51) 



Then we have 

(AA(Ax, Ay, AX U , AXi)(Z), W 3 ) = trace W 3 T AA(As)(Z) 



= ( (AA)* (-[W 3 Z + ZW[]), 



( Ax \ 

Ay 

AX V 
\AXJ 



Now we find (AL s )*(w 4 ), which consists of three columns of block with four rows per column. 
We list this by columns C\, C 2 , C 3 . 



Ci = 



/ (/C^)*(/C^ 4 )) \ 

(/c^)*(/c^K)) 

svec [{K, x Hu ){w A )\ 
\-svec [(K%)(w A )]J 



( vec (AMat (w 4 )) \ 






C 3 = (AA)*(-[Mat (w 4 )X T + XMat (w 4 ) T ]) 

Thus, the desired adjoint is given by (AL S )*(«; 4 ) = C\ + C 2 + C 3 . 

Now we evaluate (F')*(w 1 ,w 2 , W 3 , w 4 ). This consists of four columns of blocks with four 
rows per column. We list this by columns Col±, Col 2 , Col 3 , Col 4 . 



Coh 



f-()C x Hu )*(sM^(X u ow 1 ))\ 
■ -(/C^)*(sMat(A u o«; 1 )) 

Wl o s u 

V o / 
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CoU 



/(/Cy*(aMat(A,ou; a ))\ 
(/Cy*(sMat(A,o«; 2 )) 


\ w 2 os t ) 



CoU = Co/31 + Col 



'32 



CoU = (AL S 



where 



Col 



31 



/ -±vec (W 3 T AMat (x) + AW 3 Mat (2)) \ 
±svec (AW 3 + W[A) 



V 








Col 32 = (AA)*(^[W 3 Z + ZW 3 T }) 

where w x G ^ m+n \w 2 G K*( m +«), W 3 G M nxn and w 4 G K nr . Thus the desired adjoint is 
given by (F')* = Col x + Col 2 + Col 3 + ColA. 



7 Numerical Tests 

We now present results on randomly generated SNL problems with connected underlying 
graphs. The tests were done using MATLAB 7.1. The method for generating the tests 
follows from the approach used in [TTl [T9] . 

The first set of tests compares the two methods for finding a proper factorization to 
estimate the sensor locations from the optimum of the SDP relaxation. The second set of 
tests compares the two methods for solving the SDP relaxation, i.e. using the quadratic 
constraint XX T Y and the linear one using Z s >z 0. 



7.1 Two Methods for Estimating Sensor Locations 

Two methods for estimating the sensor locations from a given optimum of the SDP relaxation 
were presented in Section 14.31 i.e. 

( I X T \ 

1. Method 1: After obtaining Z s = I Jl J , X is used to estimate sensor positions 
as in [3 El HZ]. 

2. Method 2: Use the rows of P\Q to estimate the locations of the sensors. Here Q = 
VqUq, UqT,qVq = A T P 2 is the singular value decomposition for A T P 2 . 
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We denote X e i,X e2 as the estimated sensor locations from method 1 and method 2, 
respectively. We first note there is a significant difference in norm between the estimates of 
X from Method 1 and Method 2, see Table [TJ In Tables I2"f3"f4~l we use the following three 
measures to compare the two methods. 

Measure 1: Objective Function with Different Anchors 

When finding the SVD decomposition of AA T in method 2, the anchor locations are 

estimates, A e2 , and may not correspond to A. This measure uses the true objective 

Y Y T Y ■ A T 



function 



Wo[K 



A Y T A A T 



E 



1,2. 



Measure 2: Total Distance Error with True Sensor Locations 

This is a common criterion used for SNL . We compare the sum of distances between 
estimated sensor locations and true sensor locations, i.e. \\X ei — X*\\p, i = 1,2, where 
X* denotes the true sensor locations. 



Measure 3: Objective Function with Original Anchors 

We use the same criterion as in Measure 1, except that we keep the anchors fixed to 



their true locations, i.e. 



Wo }C 



X e iXTi 



AX 



T 



X ei A> 
AA T 



-E 



1,2. 





test 1 


test 2 


test 3 


test 4 


test 5 


test 6 


test 7 


mean 


std 


\\X e l — X e 2 \\f 


1.2351 


1.3002 


1.4210 


1.2906 


1.1300 


1.3810 


1.2964 


1.2935 


0.0879 



Table 1: Difference of the estimates from the two methods 





test 1 


test 2 


test 3 


test 4 


test 5 


test 6 


test 7 


mean 


std 


Method 1 


3.5200 


3.8549 


3.8793 


3.5006 


2.9434 


3.4693 


3.8736 


3.5773 


0.3363 


Method 2 


0.7462 


0.9213 


0.8656 


1.0310 


0.7237 


1.6671 


1.2351 


1.0271 


0.3319 



Table 2: Measure 1: Use X estimates in objective function; with different anchors 





test 1 


test 2 


test 3 


test 4 


test 5 


test 6 


test 7 


mean 


std 


Method 1 


1.2780 


1.4200 


1.4801 


1.3696 


1.1820 


1.4317 


1.3912 


1.3647 


0.1021 


Method 2 


0.1887 


0.1630 


0.1050 


0.1394 


0.0778 


0.0808 


0.3881 


0.1633 


0.1074 



Table 3: Measure 2: Use distance of X estimates from true sensor locations X* 



In the tests in Tables I2|3l4l we used randomly generated graphs with parameters: r = 
2, n = 16, m = 5, and radio range 0.15. The density of edges that are known was 0.75 and all 
the sensors/ anchors lie within a 2 x 2 square. We also tested many instances with different 
parameters, e.g. more sensors, and larger radio range. But the results of comparing the two 
methods were similar to those presented in these tables. 
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test 1 


test 2 


test 3 


test 4 


test 5 


test 6 


test 7 


mean 


std 


Method 1 


3.5200 


3.8549 


3.8793 


3.5006 


2.9434 


3.4693 


3.8736 


3.5773 


0.3363 


Method 2 


0.2771 


0.3264 


0.1588 


0.1799 


0.1714 


0.1453 


0.6428 


0.2717 


0.1770 



Table 4: Measure 3: Use X estimates in objective function; with original anchors 



7.2 Two Methods for Solving SNL 

In Figures fl~|2l we present results for using the quadratic constraint Y — XX T >z compared 

f I X T \ 

to the linearized version f ^ ^1^0. We solved many randomly generated problems 

with various values for the parameters. We present typical results in the figures. 

Figure [1] shows the (—log) of the optimal value at each iteration. Figure [2] shows the 
(—log) of the relative gap. Both figures illustrate the surprising result that the quadratic 
formulation is more efficient, i.e. it obtains higher accuracy with fewer iterations. This is 
surprising, since we are using a Newton based method that should be faster on functions 
that are less nonlinear. Therefore, from a numerical analysis viewpoint, it appears that the 
linear version is more ill-conditioned, as was mentioned since the constraint is not onto. In 
addition, the figures show the high accuracy that can be obtained though these problems 
are highly ill-conditioned. 

These tests provide empirical evidence for the theoretical comparison results on different 
barriers given in [T51 Hj. The results in these references show that the central path is distorted 
due to the I in the linear formulation constraint. And, the distortion increases with increasing 
dimension of the /. This agrees with our interpretation that the linear constraint is not onto, 
and the Jacobian is singular. 



8 Concluding Remarks 

In this paper, we have analyzed the well known SNL problem from a new perspective. 
By considering the set of anchors as a clique in the underlying graph, the SNL problem 
can be studied using traditional EDM theory. Our main contributions follow from this 
EDM approach: 

1. The Slater constraint qualification can fail for cliques and/or dense subgraphs in 
the underlying graph. If this happens, then we can project the feasible set of the 
SDP relaxation to the minimal cone. This projection improves the stability and can 
also reduce the size of the SDP significantly. 

In a future study we plan on identifying the appropriate dense subgraphs. (Algorithms 
for finding dense subgraphs exist in the literature, e.g. j29j [311 [25].) 

2. We provided a geometric interpretation for the method of directly using the X from 
the optimal Z s of the SDP relaxation, when estimating the sensor positions. We then 
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square for linear;triangle for quadratic 

30 i 1 1 1 1 1 — 




iteration number 

Figure 1: Comparison for two barriers; optimal value 

proposed another method of estimating the sensor positions based on a principal com- 
ponent analysis. Our numerical tests showed that the new method gave consistently 
more accurate solutions. 

3. We used the £2 norm formulation instead of the £\ norm. This is a better fit for the 
data that we used. However, the quadratic objective makes the problem more difficult 
to solve. 

In the future we plan on completing an error analysis comparing the two norms. 

4. We solved the £2 norm formulation of the SDP relaxation with a Gauss- Newton 
primal-dual interior-exterior path following method. This was a robust approach com- 
pared with the traditional symmetrization and a Newton method. We compared using 
the quadratic constraint with the linearized version used in the literature. The nu- 
merical results showed that the quadratic constraint is more stable. This agrees with 
theoretical results in the literature on the deformation of the central path based on the 
size of the I in the linearized version. 

Future work involves making the algorithm more efficient. In particular, this requires 
finding appropriate preconditioners. 
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square for linear;triangle for quadratic 




-2 1 1 1 1 1 1 1 

5 10 15 20 25 30 35 

iteration number 

Figure 2: Comparison for two barriers; relative gap 
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